In [7]:
from mpl_toolkits.basemap import Basemap
import pygrib as pb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import glob

#Get data path
datapath = glob.glob("./21072000/M-*grb2")
print(datapath)
['./21072000\\M-A0064-21072000-000.grb2', './21072000\\M-A0064-21072000-006.grb2', './21072000\\M-A0064-21072000-012.grb2', './21072000\\M-A0064-21072000-018.grb2', './21072000\\M-A0064-21072000-024.grb2', './21072000\\M-A0064-21072000-030.grb2', './21072000\\M-A0064-21072000-036.grb2', './21072000\\M-A0064-21072000-042.grb2', './21072000\\M-A0064-21072000-048.grb2', './21072000\\M-A0064-21072000-054.grb2', './21072000\\M-A0064-21072000-060.grb2', './21072000\\M-A0064-21072000-066.grb2', './21072000\\M-A0064-21072000-072.grb2', './21072000\\M-A0064-21072000-078.grb2', './21072000\\M-A0064-21072000-084.grb2']
In [43]:
"""
利用basemap套件繪製圖形,在read_grib範例中,可以得知變數的位置,這邊一樣畫2米溫度跟10米風
需要的步驟,
1.讀取grib資料,取出所需資料,包含2米溫度、10米u風、10米v風、經緯度
2.投影相關設定值(lon_0, lat_0, lat_1, lat_2)
  lon_0的key:LoVInDegrees
  lat_0的key:LaDInDegrees
  lat_1的key:Latin1InDegrees
  lat_2的key:Latin2InDegrees
3.繪製圖形
"""
import sys
cnt = 0 # 只是為了抓取1測的固定值用 
for grb in datapath:
    if cnt == 0:
        forget = pb.open(grb).select()[:][0]
        analDate = forget.analDate.strftime("%y%m%d%H")
        lon_0, lat_0 = forget["LoVInDegrees"], forget["LaDInDegrees"]
        lat_1, lat_2 = forget["Latin1InDegrees"], forget["Latin2InDegrees"]
        lats, lons = forget.latlons()
        proj = Basemap(projection="lcc",resolution='l',lat_0 = lat_0 ,lon_0 = lon_0 ,lat_1 = lat_1, \
              lat_2 = lat_2, llcrnrlat = lats[0,0], llcrnrlon = lons[0,0], urcrnrlat = lats[-1,-1], urcrnrlon = lons[-1,-1])
        cnt = 1
    basicvar = pb.open(grb).select()[:]
    fcst = str(basicvar[0].forecastTime)
    valDate = basicvar[0].validDate.strftime("%y%m%d%H")
    t2 = basicvar[62]["values"] - 273.15
    u10m_data = basicvar[66]["values"]
    v10m_data = basicvar[67]["values"]
    fig = plt.figure(figsize=(16,12))
    proj.drawcoastlines(linewidth=0.9, color='slategrey')
    ctf = proj.contourf(lons, lats, t2, cmap='coolwarm',latlon=True)
    gap = 20 #沒有用gap的話,會很慘
    proj.barbs(lons[::gap,::gap], lats[::gap,::gap], u10m_data[::gap,::gap], v10m_data[::gap,::gap] ,\
               length=5,barbcolor="navy", latlon=True)
    cbar = plt.colorbar(ctf,orientation='horizontal',pad=0.02)
    cbar.set_label("$^\circ$C",size=14)
    plt.title("Initial Time: {}".format(analDate), loc="left",fontsize=16)
    plt.title("Forecast Time:{}".format(fcst),loc="center",fontsize=16)
    plt.title("Valid Time:{}".format(valDate),loc="right",fontsize=16)
    plt.suptitle("In-Fa Typhoon\n",y=0.93,fontsize=20)
    plt.show()
    plt.clf()
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>